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Abstract 

Background: Primary osteoporosis is an age-related disease, and the main cause of this disease is the failure of 
bone homeostasis. Previous studies have shown that primary osteoporosis is associated with gene mutations. 
To explore the functional modules of the PPI (protein-protein interaction) network of differentially expressed genes 
(DEGs), and the related pathways participating in primary osteoporosis. 

Methods: The gene expression profile of primary osteoporosis GSE35956 was downloaded from the GEO (Gene 
Expression Omnibus) database and included five MSC (mesenchymal stem cell) specimens of normal osseous tissue 
and five MSC specimens of osteoporosis. The DEGs between the two types of MSC specimens were identified by 
the samr package in R language. In addition, the functions and pathways of DEGs were enriched. Then the DEGs 
were mapped to String to acquire PPI pairs and the PPI network was constructed with by these PPI pairs. 
Topological properties of the network were calculated by Network Analyzer, and modules in the network were 
screened by Cluster ONE software. Subsequently, the fronting five modules whose P-value was less than 1.0e-05 
were identified and function analysis was conducted. 

Results: A total of 797 genes were filtered as DEGs from these ten specimens of GSE35956 with 660 up-regulated 
genes and 137 down-regulated genes. Meanwhile, up-regulated DEGs were mainly enriched in functions and 
pathways related to cell cycle and DNA replication. Furthermore, there were 4,135 PPI pairs and 377 nodes in the 
PPI network. Four modules were enriched in different pathways, including cell cycle and DNA replication pathway 
in module 2. 

Conclusions: In this paper, we explored the genes and pathways involved in primary osteoporosis based on gene 
expression profiles, and the present findings have the potential to be used clinically for the future treatment of 
primary osteoporosis. 

Keywords: Primary osteoporosis, Protein-protein interaction network, Functional modules, Differentially expressed 
genes, Topological properties 



Background 

As a kind of systematic skeletal disease, osteoporosis 
results from low bone mineral density and micro- 
architectural deterioration of bone tissue which leads 
to increased bone fragility and susceptibility to fracture 
[1]. Osteoporosis can be divided into two major cat- 
egories, namely primary osteoporosis and secondary 
osteoporosis. Primary osteoporosis is a gradual process 
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and patients often cannot be diagnosed until fractures 
occur. Therefore, research on the pathogenesis of primary 
osteoporosis is urgently needed, especially regarding the 
contribution of genetic factors to its pathogenesis. 

Previous studies have found that the mutations of genes 
such as those for CSF1 (macrophage colony stimulating 
factor 1) [2] and FTH1 (ferritin, heavy polypeptide 1) 
[3], and osteoporosis-related genes such as COL1A1 
(collagen, typel, alpha 1) [4], LRP5 (low-density lipopro- 
tein receptor- related protein 5) [4], RUNX2 (runt related 
transcription factor 2) [5], WNT3A (wingless-related MMTV 
integration site 3 A) [6], DKK1 (abstract dickkopf-1) [6], 



© 2014 Li et al.; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative 
Commons Attribution License (http://creativecommons.Org/licenses/by/2.0), which permits unrestricted use, distribution, and 
reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain 
Dedication waiver (http://creativecommons.Org/publicdomain/zero/1.0/) applies to the data made available in this article, 
unless otherwise stated. 



Li et a I. European Journal of Medical Research 2014, 19:15 
http://www.eurjmedres.eom/content/1 9/1/1 5 



Page 2 of 9 



and RANKL (receptor activator of NF-kB ligand) [7] are 
related to primary osteoporosis. Low-density lipoprotein 
receptor-related protein 5 (LRP5) effects the development 
of primary osteoporosis by altering bone mineral density 
[8]. As a kind of osteoclast growth factor, CSF1 plays a 
major role in osteoporosis' development [9]. Besides, 
RUNX2 may play a key role in osteoblast cell growth and 
differentiation, and RANKL contributes to vascular calcifi- 
cation by regulating bone morphogenetic protein-2 and 
bone-related proteins [10]. Although previous studies have 
identified several potential genes and proteins as determi- 
nants of peak bone mass and susceptibility for osteopor- 
osis, the studies on primary osteoporosis are still in their 
early stages. 

In this paper, microarrays were utilized to identify differ- 
entially expressed genes (DEGs) between MSC (mesenchy- 
mal stem cell) specimens of normal osseous tissue and 
MSC specimens of osteoporosis. We used bioinformatics 
methods to construct a PPI (protein-protein interaction) 
network, and analyzed the functional modules in the net- 
work. Our work may provide a new view of osteoporosis 
and evidence for the need for further osteoporosis research. 

Methods 

All persons have given their informed consent prior to 
their inclusion in the study, and all human studies have 



been approved by the China Ethics Committee and have 
been performed in accordance with the ethical standards. 

Derivation of genetic data 

The gene expression profile of GSE35956 [11] was down- 
loaded from a public functional genomics data repository 
GEO (Gene Expression Omnibus) database which is based 
on the GPL570 [HG-U133 Plus 2] Affymetrix Human 
Genome U133 Plus 2.0 Array platform (Affymetrix, Santa 
Clara city, CA state, USA). A total of ten specimens, in- 
cluding five MSC specimens of normal osseous tissue and 
five MSC specimens of osteoporosis, were available for 
further study. 

Data preprocessing 

Firstly, we analyzed the derived genetic data by the 
Geoquery and Limma packages in R language (v.2.13.0) 
[12], and the probe-level data in CEL files were con- 
verted into probe expression matrix by RMA (robust 
multi-array average) in the Affy package [13]. Then, 
probe serial numbers from the matrix were trans- 
formed into gene names by platform R/Bioconductor 
note package of the dataset chip. Finally, as one gene 
may have many corresponding probes, the average 
value of all probe expression values was used as the 
expression value of a single gene. 




Figure 1 Dendrogram of differentially expressed genes (DEGs) by cluster analysis. The abscissa axis below represents specimen names; 
ost-op represents MSC specimens of osteoporosis; old-done represents a MSC specimen of normal osseous tissue. The abscissa axis above 
represents the clustering of specimens. The left longitudinal axis represents genes, and the right longitudinal axis represents the clustering of 
genes. Red shows up-regulated genes while green shows down-regulated genes. There are two clusters of specimen clustering: one of MSC 
specimens of osteoporosis and another of MSC specimens of normal osseous tissue; and two clusters of gene clustering: down-expression genes, 
and up-expression genes in MSC specimens of osteoporosis. 
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Table 1 Top ten enriched functions of up-regulated DEGs 



GO Terms 




Gene count 


FDR 


GO:0000278 


- mitotic cell cycle 


70 


3.51E-29 


GO:0000280 


- nuclear division 


53 


1 .68E-26 


uU.UUU/Uo/ 


- mitosis 


jo 


I .Dot-ZD 


GO:0000279 


- M phase 


63 


3.57E-26 


GO:0007049 


- cell cycle 


96 


3.81 E-26 


GO:0000087 


- M phase of mitotic cell cycle 


53 


4.30E-26 


GO:0022402 


- cell cycle process 


81 


1 .07E-25 


GO:0048285 


- organelle fission 


53 


1.35E-25 


GO:0022403 


- cell cycle phase 


67 


1 .24E-23 


GO:0051301 


- cell division 


56 


1.15E-22 



FDR: false discovery rate; GO: gene ontology. 



Screening of DEGs 

The DEGs between MSC specimens of normal osseous 
tissue and MSC specimens of osteoporosis were identi- 
fied by the samr package in R software [14]. The cut-off 
criteria were delta = 1, fold change > 2 and FDR (false 



discovery rate) was less than 0.05. To guarantee that the 
screened DEGs could represent the two types of speci- 
men, cluster analyses of DEGs were performed. 

Function and pathway enrichment of DEGs 

DAVID (Database for Annotation Visualization and Inte- 
grated Discovery) [15] online analytical tools were utilized 
to analyze the significantly enriched GO (Gene Ontology) 
terms and KEGG (Kyoto Encyclopedia of Genes and 
Genomes) pathways of up-regulated and down-regulation 
DEGs. Then the GO terms and KEGG pathways with FDR 
less than 0.05 were compared. 

PPI network construction 

The PPI pairs were acquired by directly mapping the 
DEGs to String [16] database. PPI network was con- 
structed by PPI pairs whose protein interaction scores 
were higher than 0.7. Next, the topological properties of 
the PPI network, such as degree distribution, clustering 
coefficient, average shortest path and the closeness cen- 
trality were analyzed by Network Analyzer in Cytoscape 
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Figure 2 DNA replication pathway. The red stars indicate the up-regulated genes of mesenchymal stem cells from primary osteoporosis tissues 
in this pathway. 
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software [17]. Degree distribution was computed by 
counting the number of connections between various 
proteins of the network; clustering coefficient refers to 
the average density of the node neighborhoods; average 
shortest path means the average density of shortest 
paths between all pairs of nodes; and closeness centrality 
is the sum of graph-theoretic distances from all other 
proteins in the network [18,19]. 

The protein interaction scores were calculated by the 
following formula [20]: 

S(e(x,y)) = f(diff(x),corr(x,y),diff(y)) 

= - 2 EL l0 Se(Pi) 

where diff(x) and diff(y) are differential expression assess- 
ments of gene x and gene y, respectively. Corr(x, y) repre- 
sents their correlation between gene x and gene y. Where 
k = 3, pi and p2 are the P- values of differential expression 
of two nodes, p3 is the P-value of their co-expression. 



Functional modules analysis 

Functional modules of the network were explored by 
ClusterONE in Cytoscape software [17]. Then, sub-modules 
were screened under the condition of P-value significance 
less than 0.05 and having more than ten nodes in each 
module. Finally, the fronting five sub-modules with P-values 
less than 1.0e-05 in the modularity analysis were selected 
for function enrichment using DAVID online analytical 
tools. The threshold was set as FDR less than 0.05. 

Results 

DEG analysis 

One MSC specimen of osteoporosis, which gathered in 
the same cluster as normal MSC specimens, was filtered 
and deleted after expression cluster analysis of genes. 
From the remaining four MSC specimens of osteopor- 
osis and five MSC specimens of normal osseous tissue, 
we obtained 797 DEGs, of which 660 (82.81%) were up- 
regulated genes and 137 (17.12%) were down- regulated 
genes. Clusters of DEGs are shown in Figure 1. 
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Figure 3 Cell cycle pathway. The red stars indicate the up-regulated genes of mesenchymal stem cells from primary osteoporosis tissues in 
this pathway. 
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Function and pathway annotation of DEGs 

By the analysis of GO functions and KEGG pathways of 
DEGs, we found that up-regulated genes and down- 
regulated genes were significantly enriched in different GO 
terms and KEGG pathways. Down-regulated genes were 
only enriched in function of GO:0009719 (response to en- 
dogenous stimulus), while no pathway was enriched. The 
660 up-regulated genes were significantly enriched in 50 
GO terms, and theses functional processes were mainly in- 
volved in cell cycle, material preparation for mitosis and 
DNA replication. The top ten terms are listed in Table 1. 
hsa03030 pathway (DNA replication) and hsa04110 path- 
way (cell cycle) were the only two enriched pathways of up- 
regulated genes, as shown in Figures 2 and 3, respectively. 

PPI network 

There were 4,135 PPI pairs whose reliability scores were 
higher than 0.7 and 377 nodes (47.30% of DEGs) in the 



PPI network. The network of DEGs showed a high ag- 
gregate state with two aggregate cluster sub-networks 
(Figure 4). We calculated the important topological 
parameters, such as the degree distribution, clustering 
coefficient, average shortest path and the closeness cen- 
trality of the PPI network (Figure 5), and found that the 
degree distribution of PPI network nodes followed the 
power-law of network, and had small- world network 
characteristics (short average shortest path and large 
average clustering coefficient). 

Modules analysis in the network 

High aggregation is an essential characteristic of biological 
networks, and it reflects high modularization of gene net- 
works. For distinguishing the modules which had specific 
functions and different sizes, the network was usually 
divided into relative independent sub-modules before 
analysis, and by Cluster ONE software we obtained 12 




Figure 4 Protein-protein interaction networks of differentially expressed genes. The nodes indicate genes and the edges indicate interaction 
between these genes. 
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Figure 5 Topological properties of the protein-protein interaction network. (A) The degree distribution; (B) Average clustering coefficient; 
(C) Shortest path length distribution; (D) Closeness centricity. 



models (Table 2) with more than ten nodes and P- values 
of significant less than 0.05. 

The significant enriched KEGG pathways (FDR < 0.05) 
of the fronting five modules are shown in Figure 6. Be- 
sides module 4, the other four modules all took a part in 
specified signal pathways. Module 1 participated in the 
hsa04622 pathway (RIG-I-like receptor signaling path- 
way), while module 2 participated in four pathways in- 
cluding the cell cycle and nucleotide mismatch repair 

Table 2 Modules of protein-protein interaction network 



Rank 


Nodes 


Edges 


Density 


Significance (P-value) 


1 


32 


302 


0.609 


0.000 


2 


128 


3,396 


0.418 


0.000 


3 


35 


71 


0.119 


8.949E-7 


4 


21 


47 


0.224 


1.548E-5 


5 


32 


64 


0.129 


2.427E-5 


6 


18 


19 


0.124 


3.661 E-5 


7 


22 


24 


0.104 


3.294E-4 


8 


17 


41 


0.301 


8.039E-4 


9 


12 


12 


0.182 


0.003 


10 


13 


14 


0.179 


0.004 


11 


11 


11 


0.2 


0.008 


12 


14 


50 


0.549 


0.025 



pathway. CCNA2 (cyclin A2) and CCNB1 (cyclin Bl) 
were two module 2 genes participating in the cell cycle 
pathway. Besides, the spliceosome pathway was the com- 
mon pathway enriched of modules 3 and 5. 

Discussion 

Because of the high prevalence and seriousness of osteo- 
porosis, there is an urgent need to explore the pathogenic 
mechanism of primary osteoporosis. In this study, we uti- 
lized the gene expression profile downloaded from GEO 
to explore the possible functions and signal pathways of 
DEGs and related proteins in primary osteoporosis. After 
the construction of a PPI network, sub-modules in the 
network were mined, and functions of the fronting five 
modules were enriched. 

The cell cycle pathway, nucleotide excision repair path- 
way and mismatch repair pathway were three pathways in 
which genes in module 2 were mainly involved. Cell cycle 
is the basic process of cellular activities. In the evolu- 
tionary process, the cell establishes a series of control 
mechanisms to ensure strict and orderly cell cycles that 
alternately and sequentially change. Cell cycle abnor- 
malities associated with abnormal DNA replication may 
reduce the proliferation of osteoblasts, which will lead 
to the occurrence of osteoporosis. Our work showed the 
hsa04110 pathway (cell cycle) was significantly abnormal 
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Module 1 

Pathway: RIG-l-like receptor signaling pathway 




Module 2 

Pathway: Cell cycle, DNA replication, Nucleotide excision repair, Mismatch repair 



Module 3 
Pathway: Spliceosome 




Module 5 
Pathway: Spliceosome 



Figure 6 The enriched KEGG pathways of five fronting modules in the network. Red represents nodes that did not participate in signal 
pathways in the network, while other colors imply the nodes which participated in the signal pathways. Module 2 participated in four KEGG 
pathways: cell cycle, nucleotide excision repair and mismatch repair, and DNA replication. In addition, yellow represents genes participating in the 
cell cycle pathway, while green represents genes participating in DNA replication, nucleotide excision and mismatch repair pathway as these 
pathways shared some common genes {RFC5, RFC3, RFC4, P0LE2, PCNA, and RPA3). 



in osteoporosis, which including the two down- regulated 
genes, CCNA2 and CCNB1, CCNA2 as a general regulator 
of cell cycle is the main A-type cyclin presenting in 
somatic cells and a mediator of the cell cycle [21]. This 
indicates a role for CCNA2 in the mitotic entry or com- 
pletion. Thus, the function of CCNA2 is important for 
mitosis per se, and CCNA2 could influence breakdown 
of the nuclear-envelope [22]. Moreover, CCNA2 might 
regulate some aspects of CCNB1 function [23], and 
CCNB1 is highly enriched in an ex vivo co-culture 
model of malignant breast epithelial cells and primary 
human osteoblasts by Gene Ontology analysis [24]. Pre- 
vious studies have reported that in the process of pri- 
mary osteoblast culture, cell hyperproliferation could be 
decreased and cell cycle regulation could be inhibited 
by inhibition of CCNA2 transcription [25]. So, we can 
infer that CCNA2 and CCNB1 expression represent a 
proliferative response of the osteoblast, and thus plays a 
role in primary osteoporosis. 

Our study showed that four DEGs, DDX58 (DEAD box 
polypeptide 58), IFIH1 (interferon induced with helicase C 
domain 1), ISG15 (ISG15 ubiquitin-like modifier) and 
TMEM173 (transmembrane protein 173), participated in 
RLRs (RIG-I-like receptors) signaling pathway of module 
1. DDX58 and IFIH1 are reported respectively to have four 
and two SNPs that lead to structural or polar differences 



in the mutated amino acid compared with the native resi- 
due [26]. ISG15, as the ubiquitin-like modifier, probably 
acts as one potential downstream mediator of Fzd9 (frizzled 
homolog 9) in the control of bone formation [27]. Add- 
itionally, TMEM173 is discovered to function downstream 
of mitochondrial antiviral signaling protein in the RLRs 
pathway [28], and RLRs are a family of DExD/H box 
RNA helicases which function as cytoplasmic sensors of 
pathogen-associated molecular patterns within viral RNA 
[29,30]. However, there is still yet no direct evidence 
to prove that IFIH1 and DDX58 are associated with osteo- 
porosis, and no research to support that the RLRs signaling 
pathway relates to the genesis of osteoporosis. 

We also observed that many DEGs were enriched on 
the pathway of spliceosome in modules 3 and 5, but there 
has been no literature to demonstrate the relationship 
between abnormal expression of primary osteoporosis and 
spliceosome. Ubiquity of alternative splicing in the body 
has already been proved, and most of them are involved in 
functional protein changes such as altering a proteins 
amino or carboxy terminus, or deletion of functional areas 
[31,32]. Besides, splicing variability of gene expression is a 
part of the process of cell growth and differentiation. 
However, our study showed that the spliceosome pathway 
may play some undiscovered role in primary osteoporosis 
and this warrants further studies. 
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Conclusion 

In conclusion, we analyzed DEGs profiles of primary 
osteoporosis by a computational bioinformatics approach. 
Meanwhile, we explored the DEGs enriched on pathways 
such as RLRs signaling, spliceosome, cell cycle, nucleotide 
excision repair and mismatch repair pathway. Our study 
may shed new light on the mechanism and treatment of 
primary osteoporosis by contributing a new perspective. 
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